rm(list = ls())
setwd("/Users/chengyuhang/cheng_development/r")
library(tinyarray)
library(stringr)
library(data.table)
a = geo_download("GSE10072")
group1 = ifelse(str_detect(a$pd$title,"Normal"),"normal","tumor")
table(group1)

x = geo_download("GSE40791")
group3 = ifelse(x$pd$source_name_ch1=="normal lung","normal","tumor")
table(group3)

b = geo_download("GSE31210")
setwd("/Users/chengyuhang/cheng_development/r")
getwd()
exp_b <- fread("gse31210.csv")
pdb <- fread("gse31210_pd.csv")
exp_b <- as.data.frame(exp_b)
rownames(exp_b) <- exp_b$ID_REF
exp_b <- exp_b[,-1]
exp_b <- log2(exp_b+1). #log2转换 
group2 <- ifelse(str_detect(pdb[6,],"tumors"),"tumor","normal")
table(group2)

library(hgu133a.db)
install.packages("hgu133a.db")
library(BiocManager)
BiocManager::install("hgu133a.db")
ids <- toTable(hgu133aSYMBOL)
a$exp = trans_array(a$exp,ids)

library(hgu133plus2.db)
ids <- toTable(hgu133plus2SYMBOL)
exp_b = trans_array(exp_b,ids)
x$exp =trans_array(x$exp,ids)

y = read.table("TCGA-LUAD.htseq_fpkm.tsv.gz",
               header =  T,
               check.names = F,
               row.names = 1)
group4 = make_tcga_group(y)
table(group4)
y = trans_exp(y)
y = y[rowSums(y)>0,]

genes = intersect_all(rownames(a$exp),
                      rownames(exp_b),
                      rownames(x$exp),
                      rownames(y))
length(genes)

install.packages('VennDiagram')
library("VennDiagram") #看重复的RNA
draw_venn(x = list(a = rownames(a$exp),
                   b = rownames(exp_b),
                   x = rownames(x$exp),
                   y = rownames(y)),
          "genes")

library(AnnoProbe) #看看是什么类型的RNA
anno = annoGene(genes,"SYMBOL")
table(anno$biotypes)

exp = cbind(a$exp[genes,],
            exp_b[genes,],
            x$exp[genes,],
            y[genes,])
exp = as.matrix(exp) #转换为矩阵更好进行运算 而不是数据集
dim(exp)

Group = c(group1,group2,group3,as.character(group4))
Group = factor(Group,levels = c("normal","tumor"))
table(Group)

mod = model.matrix(~as.factor(Group))
batch = rep(1:4,time = c(length(group1),length(group2),
                         length(group3),length(group4)))
table(batch)

library(sva)
install.packages("sva")
BiocManager::install("sva")
exp_adj = ComBat(exp,batch = batch,mod = mod,par.prior = TRUE,ref.batch = 1)

fids = data.frame(probe_id = rownames(exp),
                  symbol = rownames(exp))
degs <- get_deg_all(
  exp,
  Group,
  ids,
  logFC_cutoff = 2,
  pvalue_cutoff = 0.05,
  adjust = TRUE,
  entriz = FALSE,
  
  n_cutoff = 3,
  cluster_cols = TRUE,
  annotation_legend = FALSE,
  show_rownames = FALSE,
  
  lab = NA,
  
  symmetry = FALSE,
  
  color_volcano = c("#2874c5","grey","#f87669")
)
degs$plots


